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ABSTRACT 



We present three dimensional simulations of the interaction of a light hypersonic jet with an 
inhomogeneous thermal and turbulently supported disk in an elliptical galaxy. These simulations 
are applicable to the GPS/CSS phase of some extragalactic radio sources. 

We identify four generic phases in the evolution of such a jet with the interstellar medium. The 
first is a 'flood and channel" phase, dominated by complex jet interactions with the dense cloudy 
medium close to the nucleus. This is characterized by high pressure jet gas finding changing weak 
points in the ISM and flowing through channels that form and reform over time. A spherical, 
energy driven, bubble phase ensues, wherein the bubble is larger than the disk scale, but the 
jet remains fully disrupted close to the nucleus, so that the jet flux is thermalised and generates 
a smooth isotropic energy-driven bubble. In the subsequent, rapid, jet break-out phase the jet 
breaks free of the last obstructing dense clouds, becomes collimated and pierces the more or 
less spherical bubble. In the final classical phase, the jet propagates in a momentum-dominated 
fashion similar to jets in single component hot haloes, leading to the classical jet - cocoon - 
bow-shock structure. 

Subject headings: Radio Galaxies — GPS/CSS ISM — Turbulence, Fractal medium. 



1. Introduction 

In this paper we present a high resolution simulation of a radio jet interacting with an inhomogeneous 
interstellar medium in the form of a turbulently supported disk which extends our previous work in several 
different directions: (1) The simulation is three dimensional with a resolution of 512 cells in the (cartesian) 
coordinate directions and a spatial scale of 2 parsecs per cell. (2) The density structure of the warm disk 
is described by a log-normal distribution, typical of the warm interstellar medium in a number of different 
environments. (3) The initial distribution of the warm medium is that of a thick almost Keplerian disk 
supported by a combination of thermal pressure and a, dominant, supersonic velocity dispersion. This type 
of distribution is one of several that may be contemplated and is supported by the observation of such disks 
in M87 flDopita et Z1|1997[ ) and NGC 7052 ( |van der Marel and van den Bosch||1998[ ). (4) The soft X-ray 



emissivity and surface brightness of the thermal plasma is calculated, in conjunction with the radio surface 
brightness image. These synthetic images provide valuable insights for the interpretation of observational 
data. 

This simulation illustrates four important phases in the evolution of a young radio galaxy: (1) An initial 
"flood and channel" phase wherein the radio source is beginning to be established, but is still interacting 
strongly with the disk material. (2) A quasi-spherical jet-driven bubble phase where the jet is fully disrupted. 
(3) a rapid jet establishment phase where the last obstruction is ablated away and the jet reforms and crosses 
the spherical bubble. (4) A classical jet - cocoon - bow-shock phase familiar from previous studies of radio 
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sources in a single component, hot ISM. We also determine the X-ray morphology associated with the 
interaction of the radio source with the disk and show that the luminosity from the disk may persist well 
after the initial flood and channel phase . 



Model Parameters 



2.1. Host galaxy potential 



Elliptical galaxies contain both baryonic and dark matter. The former dominates on scales close to the 
core and the latter dominates at radii ~ lOkpc. Since jet simulations frequently extend over such scales, we 
have constructed a family of potentials which we use to represent a combined, self-consistent distribution 
of baryonic and dark matter. In the simulation presented here, the scale of the simulation is such that the 
potential is dominated by the luminous component close to the core. The potential is based on isothermal 
distributions of baryonic and dark matter, described in terms of isotropic distribution functions. 

Let (Td and t\d be the velocity dispersion and core radius of the dark matter distribution. Normalising 
the potential by <r 2 D and the radius by t\d, the dimensionless version of Poisson's equation is: 
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It is convenient to take the ratio of dark to luminous velocity dispersions, k, and the the ratio of core radii, 
A, as defining parameters of the potential and mass distribution. The central densities are given by: 
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We numerically integrate ([lj for chosen parameters, tabulate the results and interpolate when required. In 
this paper we use k = 2, and A = 10, referring to the potential as ^2,10 ■ The central density is dominated 
by baryonic material over dark matter by a ratio of 25 : 1. 

We adopted this isothermal formalism for the galactic potential for two reasons. First, it seems to be 



supported by observations of at least some classes of galaxies, ( c.f. Cowsik et al. (20071, Kuzio de Naray 



et al. (20061 , de Blok (2005 1 ) although perhaps not all (e.g. Schmidt and Allen (2007) ). Secondly, the well 



known Navarro, Frenk & White, (NFW), (Navarro, Frenk, and White 19971, profile arises in an empirical 
manner from smooth particle cosmological simulations, whereas the isothermal formalism arises from simple 
physical considerations (King 19661. In the absence of a definitive proof cither way, we choose to follow 



the isothermal path. We note that with a coupled baryonic/dark matter potential, the region between the 
baryonic core radius (350 pc) and the outer dark matter core radius (3.5 kpc) does produce intermediate 
density radial power-law indices in the range; < a < — 1, before turning down to steeper slopes beyond the 
dark matter core radius (well outside our present simulation) . If a galaxy with such a double potential were 
observed with insufficient resolution to resolve the inner 300 pc, then the observed distribution inside the 
dark matter core radius would be observed to be steeper than a single isothermal distribution, perhaps being 



consistent with a more 'cuspy' NFW distribution. (See Sutherland and Bickncll (2007) for more details) 
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2.2. The Hot Galactic Atmosphere 

The tenuous, hot interstellar medium in our simulations is isothermal. For gas temperature T the dark 
matter virial temperature is T» = ma^/k, where fh = pm u « 1.035 x 10 -24 g, is the mean particle mass. 
Assuming hydrostatic equilibrium, the particle number density, of the hot gas is nh/«h,o — exp [—T*/Tip]. 



2.3. The non-uniform, fractal, warm interstellar medium 



To establish a non-uniform interstellar medium we utilise an obvious analogy with a turbulent medium. 
Rather than attempt to review the vast topic of turbulence here, we refer the reader to the recent astrophys- 
ically oriented reviews of (Elmegreen and Scalo 2004), and ( Scalo and Elmegreen||2004 1 as a starting point 
and then refer simply to some specific results that we use below to arrive at reasonable parameters for our 
non-uniform medium. 



2.3.1. Log-normal density distribution 

We use a log-normal distribution to describe the single point statistics of the density field of our nonuni- 
form ISM. The log-normal distribution is a skewed continuous probability distribution, for an independent 
variable x ^ 0. Unlike the normal distribution, it has a non-zero skewness, variable kurtosis/flatness, and 
in general the mode < median < mean. The log-normal distribution appears to be a nearly universal prop- 



erty of isothermal turbulent media, both experimentally, numerically and analytically ( e.g. Nordlund and 



Padoan (19991 , see also |Warhaft (2000) ). It is also at least intuitively encouraging that is corresponds 



the to the limiting distribution for the product of multiplicative random increments, in the same way that 
the normal distribution plays that role for additive random increments in the central limit theorem. It is 
thus compatible at least conceptually with a generic cascading process consisting of repeated folding and 
stretching. 

The probability density function with parameters s and m (the standard deviation and mean of the 
corresponding normal distribution), for the log-normal ISM density p is, 

P(p) = l/(sV27rp)exp [-(lnp - m)/2s 2 ] . 



We adopt a mean 



and 



fi = exp[m+ s 2 /2] = 1.0. 



a 2 = p 2 (exp[s 2 ] - 1) = 5.0, 



as our standard log -normal distribution, compatible with the favored ranges in Fischera, Dopita, and Suther 



land (20031 and Fischera and Dopita (20041 in their star burst galaxy reddening/exctinction models. The 



variance determines the concentration of mass is into dense cores, or conversely the volume occupied by 
voids. With these parameters, densities below the mean, p, comprise one quarter of the mass, and three 
quarters of the volume; the mean is approximately 20 times the mode. 
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2.3.2. Powerlaw density structure 



The isotropic power spectrum of the fourier transform of the density, F(k) is D(k) = J k 2 F(k)F* (k)dH,, 
where the angular integrals are over the relevant angular variables in Fourier space. If D(k) oc fc _/3 , and 
if P — 5/3, the spectrum is often referred to as Kolmogorov turbulence. A scalar tracer (density) of the 
turbulent field also shows the Kolmogorov structure index (War haft 20001; hence we generate a density 
distribution with the Kolmogorov index for our inhomogeneous ISM. 



2.3.3. Equilibrium turbulent disks 



Many radio galaxies exhibit a turbulent disk of gas in the central regions. We establish equilibrium 
turbulent disks in the appropriate potential by first determining the ensemble average distribution of the gas 
and then realizing a single instance of the ensemble by multiplying by a log-normal power-law structure as 
described above with unit mean, variance a 2 — 5.0 and power-law index (3 — 5/3. 



Following, Strickland and Stevens (2000) we adopt the ansatz of an almost Keplerian disk. We find 

e 2 Vv,o - (1 ~ e 2 )^ ,o} 
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where o~ g is the turbulent velocity dispersion of the gas. The main differences from the development of a 



similar equation by Strickland and Stevens (2000) is the (logically required) lack of dependence of e# on 



z and the formal introduction of a turbulent velocity. The latter avoids the difficulty of prescribing an 
unphysically large temperature in order to achieve a reasonable disk scale height. 

In summary we use the double isothermal potential ip K ^\with k — 2, A = 10, scaled by ru = 3.5 kpc, 
and ajy = 400 km/s, and a rotation parameter ex = 0.93, plus a velocity dispersion of the warm gas of 
dg = 40 km/s for the disk. 



Jet Parameters 



We use a non-relativistic code for these simulations. The use of such a code to simulate phenomena 
that involve relativistic flow in parts of the grid, is not ideal. Nevertheless, a large part of the flow field is 
non-relativistic and is driven by the energy or momentum flux provided by the jet. Hence, we establish jet 
parameters in such a way that the energy flux of the non-relativistic jet corresponds to the jet energy flux 



of a given relativistic jet using relationships derived by Komissarov and Falle (1996) 



The density ratio, r\ and the Mach number, M nr of a non-relativistic jet corresponding to a jet with 
Lorentz factor T and density parameter \ = pc 2 /4p is described by the following equations. The parameters £ 
and T; sm are the ratio of jet to external pressures and the (hot) interstellar medium temperature respectively. 
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(4) 
(5) 



Note that the low value of kT- lsm /mc 2 ~ 10 6 guarantees a light non-relativistic jet (i) < 1) and that the 
non-relativistic Mach number effectively corresponds to the Lorentz factor. 
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In our simulations we take the gas to have the ideal adiabatic index, 7 = 5/3, since this represents the 
external medium the most accurately and we are mainly interested in the effect that the jets have on the 
external clumpy medium. 



3.1. Resolution and Scaling 

We have selected a jet for which the instantaneous hot spot /3hs ~ 0.05 — 0.1 and the initial jet diameter 
is 20 pc. These parameters, together with a central ISM p/k — 10 6 cm -3 K places the jet at the low end of 
FR2 power « 3 x 10 43 ergs s _1 . This choice of parameters produces a jet which initially interacts strongly 
with the ambient ISM but whose morphology at later times is similar to a classical FR2 radio source. 

With adiabatic simulations, the choice of the spatial, velocity and density scales is arbitrary. However, 
the introduction of cooling (in the thermal gas) restricts the allowable scaling to a one-parameter set. Here 
a given simulation defines a one-parameter family of simulations where the scaling parameters satisfy the 
following constraints: 

xq = arbitrary, vq — fixed , 

t = xq/vq, p cx x„ 1 , ,g. 

Po = PqvI, kT /m = v%, 

p*,o oc Xq 2 , M* >0 oc x Q . 

with the proviso that xq cannot increase to the extent that the gas density exceeds the density of gravitating 
matter. Notwithstanding this restricted one-parameter scaling, the allowable set of physical models allowed 
by this scaling describes an interesting variety of different physical situations. (See Sutherland & Bicknell 
2007 for further details) 

Three simulations were performed, models A, B, and C. The 'Standard Model' A) is presented in the 
greatest detail. Models B and C are used to compare the effect of the density and the distribution (smooth 



or fractal) of the warm gas. Table 3.1 summarizes all the parameters used in Model A, for the jet, the 



potential and the interstellar medium respectively. The fiducial scaling uses a spatial scale of xq — 1 kpc. 



4. Slice plots and surface brightness images 

We now present multi-panel snapshots from significant epochs; these snapshots are designed to bring 
out the relevant physics of the simulations. Some snapshots are slices of important variables such as density 
and pressure; in some cases we also present projected versions of variables such as the density. We also 
present images of the radio and X-ray surface brightness obtained by ray-tracing of the output volumes. 

In all cases the times chosen for the sequence of snapshots corresponds to the following phases in the 
simulation. The mid-plane density slices (Figure [T]) most clearly illustrate the phases enumerated here, with 
the other representations highlighting some specific aspects. 

1. Flood and Channel phase: Snapshots at 5 - 15 kyr represent the time over which the jet makes its way 
through the porous, fractal disk. The shape of the interacting region is amorphous, and determined 
by flow of hot high pressure gas along weaker line in the dense disk medium. The pressure in this 
phase is very high, and the x-ray emission is a strong function of time, depending on a combination of 
the amount of disk material that is advected and the amount of energy that is processed by radiative 
shocks, which increases with the size of the region. 
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Parameter & units Value 
Equivalent Relativistic Jet Parameters 

tLorentz factor 5 

tRcst energy density /enthalpy 10 

Velocity / Speed of light 0.9798 
Hydrodynamic Jet Parameters 

Pressure / External pressure 1.0 

Density / External density 2.0 X 10 -3 

Mach number 25.9 

tDiameter (pc) 40 

Kinetic luminosity 2.77 X 10 43 



Double isothermal potential 

tDark matter to 

Baryonic velocity dispersion 
tDark matter to 

Baryonic core radius 

Central Baryonic density to 

Dark matter density 
Dark matter values 
tVclocity dispersion (kms -1 ) 
tCorc radius (kpc) 

Central density (gem -3 ) 
Enclosed Masses, at r = 

Dark Mass (Mq) 

Baryonic Mass (Mq) 

Total Mass (Mq) 

Baryonic/Dark Mass ratio 
Baryonic values 

Velocity dispersion (kms -1 ) 

Core radius (kpc) 

Central density (gem -3 ) 
Enclosed Masses, at r = rj, 

Dark Mass (Mq) 

Baryonic Mass (Mq) 

Total Mass (Mq) 

Baryonic/Dark Mass ratio 
Hot Atmosphere: 
tVirial/Gas temperature 
Gas Temperature (°K) 
Central Values: 
tpressure/fc (cm" 3 °K) 

pressure (dynes cm -2 ) 

number density (cm -3 ) 

mass density (g cm -3 ) 
Warm Disk-ISM : 



Virial/Gas temperature 1200.0 

fGas Temperaturc(°K) 1.0 X 10 4 

tTurbulcnt dispersion (kms -1 ) 40.0 

tRotational Support 0.93 
Internal Non— Uniformity : 

tLog-Normal Mean 1.0 

fLog-Normal Variance 5.0 

fDcnsity Power Law 5/3 

Volume of warm gas (pc 3 ) 2.55 X 10 7 

Mass of warm gas (Mq) 4.67 X 10 5 
Central Values: 

pressure/fc (cm -3 °K) 1.00 X 10 6 

tnumbcr density (cm -3 ) 10.0 

mass density (g cm -3 ) 1.04 X 10 -23 



Assigned parameters are indicated with a tsymbol 

Table 1: Standard Jet, Host and ISM Parameters. 



10 

25 

400 
3.5 

1.47 x 10 -2; 

8.14 x 10 10 
4.56 x 10 10 
1.27 x 10 11 
0.6 

200 
0.35 
3.68 x 10 -2; 

3.08 x 10 s 
4.72 x 10 9 
5.03 x 10 9 
15.3 



1.0 
1.20 x 10 7 

1.00 x 10 6 
1.38 x 10 -10 
8.35 x 10 -2 
8.64 x 10 -26 
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2. Energy Bubble phase: Snapshots at 25 - 45 kyr represent the epoch during a high pressure bubble 
has grown larger than the disk, to form a smooth nearly spherical bubble in the hot atmosphere. The 
jet is still disrupted in the disk, but the bulk of the energy flux drives the expansion of the nearly 
adiabatic bubble; there is a corresponding drop in the efficiency of conversion of jet energy to thermal 
emission. Most of the bubble grows (outside the disk) with a Radius-time power law that is consistent 
with classical energy driven bubble theory. The flood and channel behavior continues in the disk plane, 
but this involves a steadily decreasing fraction of the jet energy flux. 

3. Jet Breakout phase: In the 55 - 70 kyr epoch the jet breaks free of the few remaining clouds in its 
path and the jet terminus propagates towards the edge of the bubble. The Bubble has a low density, 
so that the jet transits the bubble quickly. 

4. Classical Phase: At t 75kyr the jet pierces the original bubble and subsequently forms a classical radio 
lobe with hotspot, cocoon, bow-shock and backflow. The remnant of the spherical bubble continues to 
grow, and flow continues within the disk, with numerous radiative and non-radiative shocks throughout 
the whole disk persisting to late times excepting the very centre of the disk, which has been cleared 
by the main jet. 



4.1. Mid-plane Density and Pressure Slices 

Figure [l] shows a sequence of density slices. At 10 and 15 kyr "bright" spots of high density, caused 
by radiative shocks driven into the warm gas, are evident. The jet is anything but straight as it selects a 
path of least resistance through the region in which the initial density is described in terms of a log-normal 
distribution (see § s:lognormal). 

In the (25 - 45 kyr) phase a dense obstruction causes the jet to split. Further regions of high density 
are observed, often at the tips of dense cloud material that is being ablated by the outflowing radio plasma. 
The thermalized jet pressure drives an almost spherical bubble in the surrounding medium. A bow-shock 
surrounds the bubble and the contact discontinuity between the bubble and the hot interstellar medium is 
apparent. 

In the jet-breakout phase (55 - 70 kyr) the pressure of shocked jet plasma, continues to drive a more 
nearly spherical bubble into the hot interstellar medium. Towards the end of this phase, the obstructing 
cloud that was responsible for the disruption of the jet dissipates and a straight jet channel is starting to be 
established. 

At 75 kyr (the lowest left panel of Figure [T]) a unidirectional jet flow has been established and the jet is 
about to pierce the initial bubble. Subsequently, the jet propagates beyond the bubble and starts to establish 
a classic FR2 lobe morphology consisting of bow-shock, cocoon and backflow. Another way of looking at 
this phase is that the quasi-isotropic (energy-driven) phase of radio lobe expansion makes a transition to 
a bow-shock dominated (momentum-driven) phase as the now relatively undeflected jet progresses through 
the interstellar medium. An interesting feature, which is discussed further below, is that the jet, whilst 
maintaining a more or less single direction, is unstable. As a result of both filamentation and helical 
instabilities, the end of the jet has become somewhat diffuse and the jet momentum is spread over a region 



which is wider than the original jet diameter in a similar manner to that envisaged by Scheuer ( 1982 1 



Mid- plane pressure slices in Figure [2] confirm the evolutionary features revealed by the density and at 
the same time show some additional characteristics. 
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Fig. 1. — Mid-plane density slices: The panels represent the logarithm of the density in a k = 255 plane of the 512 3 simulation 
in the various phases of the evolution described. The grayscale bar shows the range of the density. The selection of snapshot 
times is the same in all subsequent figures. 
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The disk interaction phase shows a high-pressured amorphous region where the jet struggles to make 
its way out of the confining disk. The lower pressure dense gas (the white region) surrounds this amorphous 
bubble. In the flood and channel phase the pressure in the jet-driven bubble is fairly uniform (as expected) 
but also shows at 35 kyr a biconical shock associated with the compression of the emerging jet by the 
overpressured bubble together with a bow shock caused by obstructing high density gas in the path of the 
jet. There are also some spots of high pressure associated with the general turbulence in the bubble. During 
the jet breakout phase most of the radio lobe is at constant pressure, and there are only small residual 
regions of low pressure associated with the disk. Further biconical shocks form within the jet and the high 
pressure jet terminus is seen propagating through the bubble. During the classical phase the hot spot at 
the end of the jet pierces the bubble and forms a higher pressure, momentum-driven cocoon. The biconical 
shocks do not persist within this cocoon and the high pressure region at the end of the jet is spread out over 
more than a jet diameter. 

4.2. Radio Surface Brightness 

We construct a surface brightness emissivity from the pressure p n t of non-thermal plasma, defined to 
be that which originates from the jet. The emissivity is defined by: j v oc v =a where p n t is the 

nonthermal pressure and a is the spectral index (with the convention F v oc v~ a ). This expression assumes 
that the magnetic pressure tracks the plasma pressure, that is, B 2 /8tt oc p. The scalar tracer phi, which is 
the mass concentration by density of the non-thermal plasma, is used to indicate the presence of jet plasma 
at the various regions of the grid. We use the above expression for the emissivity when <f> > 1CP 3 . 

A series of synthetic radio surface brightness snapshots is displayed in Figures [3] and [3] at the same times 
as in the previous figures. At the early times the surface brightness shows the amorphous structure that is 
evident in the density slices discussed previously. 

In the flood and channel phase two radio-emitting bubbles form. There is some flocculent structure in 
the surface brightness resulting from the non-uniform thermalisation of jet plasma during this phase. We 
also see quite clearly (particularly at 35 kyr and 45 kyr) the bow shock (on each side) caused by the impact 
of the jet on the obstructing cloud discussed in § |4.1| 

During the jet-breakout phase the radio-emitting bubbles become smoother in appearance but high 
surface brightness features associated initially with the jet-cloud bow shock (55 kyr) and then with the jet 
terminal shock (65 and 70 kyr) are also seen propagating out through the respective bubbles. Another 
feature which becomes obvious here, but which is also evident from approximately 25 kyr onwards, is a band 
of low surface brightness emission caused by the partial exclusion of radio-emitting plasma from the disk 
region. 

From 75 kyr onwards we see the jets breaking free of the bubbles. 

4.3. X-ray emission 

Images of the X-ray surface brightness present a different view of the physics of this simulation. The 
temperature of this gas and the total thermal cooling per unit volume, p 2 A p (T), determine the spectrum of 
gas at each voxel, via the original MAPPINGS spectrum used to define the cooling function. The spectrum 
is integrated over specific bands to determine X-ray emissivities in those bands. These emissivities are then 
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Fig. 2. — Mid-plane pressure slices: The panels represent the logarithm of the pressure in an k = 255 plane of the 512 3 
simulation at the same phases of the evolution as in Figure [l] Note that for clarity, the sense of the grayscale colorbar has been 
reversed compared to that for the density. 
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Fig. 3. — Radio Surface Brightness: The panels represent the logarithm of the synthetic radio surface brightness at the 5 — 
45 kyr phases of the evolution, corresponding to the first 6 panels of the density evolution in Figure [l] The insets for 5 kyr and 
10 kyr show the surface brightness images magnified by factors of 3 and 2 respectively in order to bring out the detail in these 
early stages. 




Fig. 4. — Radio Surface Brightness: The panels represent the logarithm of the synthetic radio surface brightness at the 55 — 
95 kyr phases of the evolution, as in the second 6 panels of Figure^ 
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integrated along the line of sight to determine the X-ray surface brightnesses. With three bands, false colour 
RGB images can be constructed and these are presented in Figures [5j and [6] 

The soft X-ray image, Figures [5] seen nearly face on to the disk at t = 50 kyr, uses bands between 0.1 
and 1.5 keV for red and green components, with the hard 1.5-10keV component in (faint) blue. This view 
shows how the dense clouds show sub-structure in the X-rays, where red dots pick out slow fully radiative 
Shockwaves at the limit of the simulation resolution. The higher energy green and blue cmisison traces the 
hotter non-raditatively shocked ISM gas that is subsequently advected along with the hot bubble. This 
radiative/non-radiative divide separates inefficient and efficient mass transport, and a significant fraction of 
the densest gas is radiatively shocked and only inefficiently transported by the expanding bubble. Figures [6] 
suggests that the majority of the advected gas and the hot bubble itself is primarily visible in in emission 
above 2 keV. 

The total 0.1 — l.OkeV X-ray emission for models A, B and C as a function of time is shown in Figure|7] 
The common feature of all the luminosity curves is that they rise approximately linearly to a broad maximum 
and then decline. For the more realistic models A and B, the maximum luminosity is a fraction of one percent 
of the jet energy flux. In the physically unrealistic model C, the peak X-ray luminosity is about 1.6% of 
the jet energy flux. In all cases, the peak X-ray luminosity occurs when the jet can be considered to have 
just broken free of the disk. Before that epoch a portion of the jet power is being directed into the radiative 
shocks that are responsible for the X-ray emission. When the jet breaks out, the fraction of its power that 
was being diverted into radiative shocks is reduced since the jet power is now diverted into non-radiative 
shocks in the hot interstellar medium. 



5. Mass Redistribution 

One of the main motivations for examining jet-ISM interactions in detail is to examine the importance 
of black hole induced feedback processes in galaxy formation. Here we consider the mass that is driven to 
larger scales as a result of the interaction of the jet with the turbulent disk. 

In addition to the main simulation model A, models B and model C were performed to look for changes 
in the interaction sequence with a change in two disk ISM parameters, density and uniformity. Model B 
differs only in the mean density of warm gas. The central value in model B is 20 cm -3 compared to 10 cm -3 
in model A. The point of model B is to examine the effect of the mean density of identically distributed gas. 
Model C, has the same mean density as model B, but has a perfectly smooth distribution of warm gas. The 
point of this model is to examine the effect of the porosity of the gas distribution. 

The resolution of Models B and C is a factor of two lower, i.e. 256 x 256 x 256. However, this does 
not appear to affect the comparison; similar behaviour is observed and the jet is adequately resolved with 
10 resolution elements across its diameter. We do not expect the much more uniform density distribution in 
model C to suffer from lower resolution. 

We have divided the computational domain into a nested series of rectangular sub-zones, and define 
integration regions as the difference between successive zones as follows. Integration of the density in each 
zone give the zone masses as a function of time, and the differences between successive zone integrals allow 
us to monitor the transport of material from region to region. The zones are rectangular for integration 
simplicity, although the regions between them may not be. Figure [8] shows the regions graphically for a late 
time density snapshot of Model A. 
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Fig. 5. — Soft X-rays: the logarithm of the soft X-ray surface brightness nearly face on at 50 kyr. Red indicates emission 
between 0.1 and 0.75 keV, and marks the location of slow, fully radiative shocks. Green represents emission between 0.75 and 
1.5 keV, and comes from gas > 10 7 K, while blue comes from emission betwen 1.5 and 10 keV, and shows how weakly the very 
hottest gas (> 10 9 K) is cooling. 




Fig. 6. — Hard X-rays: as for the previous figure, but for the hard X-ray bands Red = 2—4, Green = 4-6 and Blue = 6-10!keV 
emission. The very hot outer bubble shock is just visible in the hardest bands, and the red emission shows gas up to ~ 3 X 10 7 K 
being advectcd around the dense disk clouds. 
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Fig. 7. — The evolution of the 0.1 - 10.0 keV X -ray luminosity of models A , B and C. The times of the snapshots used in 
the detailed description of model A are indicated by filled circles. A line indicating 1% of the jet energy flux is also drawn. 
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Mass redistribution in the simulations is summarised in Figures [9] and 10 The upper panels in Figure [9] 



show, for each model, the disk mass fraction, defined as the fraction of mass in each region occupied by 
matter originally in the warm disk, as a function of time. The lower panels show the time rate of change 
of mass in each region as a function of time. Region 3 shows a range of processes, uplift from zone 1 and 
fallback to zones 1 and 2, and the mass in region 3 is variable. The upper atmosphere regions, 5 and 6, 
see very little of the disk material by the end of the simulations, and even with significant movement of 
the disk in model bf C, 2% or less ends up above 500 pc from the disk plane. The mass exchange rates in 
models A and B are similar, peaking at around 2 M Q per year in the plane of the disk between region 1 
and 2, while the monolithic disk shock in the uniform test model C generates rates up to 5 times greater. 
Models B and C are identical except for the non-uniformity of B suggesting that quantitative models of mass 
transfer in jet-host feedback models needs to take non-uniformity into account. Detailed knowledge of the 
non-uniformity in real radio galaxy hosts is needed before truly realistic transport models can be computed. 

Figure [10] shows a projection of the density in the first 3 regions at selected epochs. It illustrates how 
the single disk shock in the uniform model, C, is effectively sweeping the inner regions of the disk clear of 
high density gas, in contrast to the non-uniform models A and B. 

Three key points about the mass redistribution are evident from the figures: 

1. Mass Advection Efficiency with a Radiative Critical Density The initial and final mass frac- 
tions in Figure [9j show that the two clumpy models (A and B) are similar, with a significant fraction 
of disk material remaining in the inner region 1. The relative mass remaining in Model B is some- 
what greater and A. This may be consistent with the notion of a critical density, p CI , above which the 
traversing shocks become locally fully radiative, cool and compress gas efficiently forming small dense, 
slow -moving clumps. For p < p CY i t the shocks may be adiabatic, heating and expanding the gas which 
becomes advected rapidly along with the hot gas of the main blast. 

If there is such a critical density p cr , then we expect that in the denser Model B the mass fraction 
compressed into dense clumps would be higher since p cr is determined by the jet ram pressure and the 
cooling function, and the corresponding fraction above a fixed point in the denser model is greater. In 
region 1 in Model A, about 60% of the original mass remains at late times, and in model B at a similar 
late dynamical time (estimated by the extent of the main jet-bubble), 72% remains. This suggests that 
as long as the material above a radiative critical density remains in dense, un-connected, clumps, the 
non-gaussian log-normal statistics of the density may influence the behavior of the global system. 

2. Inefficient Transport of Disk ISM to outer Hot Atmosphere In all cases, about 83-85% of 
the gas starts in zones 1 and 2 combined, and 75-78% remains in the two zones at the end of the 
simulation, with only about 4-7% of the disk gas being transported to the higher regions 4, 5 and 
6. Only trace amounts of disk material are transported to region 6 and this is largely in the form of 
material diffused into the very hot radio plasma. This fraction is uncertain since it may be dominated 
by numerical diffusion rather than physical transport. Despite the fact that the main blastwave crosses 
the entire disk, and locally many shock crossing timcscalcs pass, the majority of the warm clumpy disk 
ISM is not moved out of the plane of the disk, and the jet outburst is inefficient at 'clearing out' the 
disk. 

3. Transport of Disk ISM within the Disk In the uniform model C, 75% of the disk mass initially 
in region 1 is swept to region 2, and as the circular empty region increases we expect that essentially 
all of the inner disk would be swept clear in another 20 — 40 kyr. In models A and B only 30-40% of 
the disk in region 1 is transported to region 2, and the timescale for clearing out region 1 is becoming 
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Fig. 8. — A 3D perspective and side view of the arrangement of the regions in the mass integration analysis. Each zone 
encompasses the preceding zones, each being a rectangular integration box, and the labeled regions are defined at the difference 
between each successive pair of zones, with region 1 and zone 1 being the same. 
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Fig. 9. — Mass fractions and derivatives in regions 1-6, for model A (left), B (center) and C (right). The upper panels show 
the fraction of disk material, by mass ,in each of the 6 regions, labelled by region number. The letters U, P and F indicate the 
dominant direction of mass transfer on the curves. U indicate periods of "uplift' . P indicate in plane movement. F indicates 
where material is predominately falling back to a lower zone. 
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Fig. 10. — These images show the distribution of gas density at the beginning, midpoint and final state for the three models 
in zone 4, i.e. just the lower 256 pc (1/4) of the grid, excluding the upper atmosphere regions 5 and 6 for clarity. The endpoints 
were determined by halting the simulation once the jet reached the far edge of the computational grid. The location of regions 
1-3 are shown. 
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greater than 150,000 years (assuming constancy of the final mass loss rates). It is therefore likely that 
dense material will remain in region 1 for a longer time in models A and B compared to model C. 

The clumpy nature of the models A and B significantly extend the lifetime of material in the inner zones 
of the outbursts, showing that in the cirumnuclear region of an active galactic nucleus, the inhomogenity of 
gas in the inner region makes it possible to retain gas there, maintaining a supply of accreting gas to the 
black hole for longer than models with homogeneous gas distributions may suggest. The relevant timescale 
for the duration of the outburst is determined primarily by the accretion timescale. 

6. Summary 

The main simulation that we have presented, of a jet interacting with a turbulcntly supported disk, 
exhibits four interesting phases: 1) A 'flood and channel' phase 2) A nearly adiabatic (globally) energy 
driven bubble phase 3) A jet breakout phase 4) A classical phase. These are described in detail in § [4] 

The late-time morphology of the radio source manifests a clear signature of the two earlier phases. 
The X-ray emission also presents some further insights. Not only do we see X-ray emission from the disk 
associated with radiative shocks in the early stages and X-ray emission associated with the bow-shock of the 
radio lobes but we also see persistent X-ray emission from the disk well after the radio bubble has passed 
through the disk. This is caused by the continued driving of radiative shocks into the disk by the high 
pressured bubble. The peak in the X-ray luminosity occurs during the jet breakout phase, when more of the 
jet energy flux is diverted into the growing, hot adiabatic bubble and there is less available for processing 
in the dense disk material. The behaviour of the disk is suggestive of a critical phenomenon: The critical 
cloud density at which the ram pressure of the jet is converted into fully radiative shocks, and the mass 
fraction that is available for rapid mass advection may be influenced by the log-normal density field in the 
disk. Additionally , the global movement of disk material throughout the simulation is clearly very different 
for a uniform disk model compared otherwise similar, albeit somewhat ad hoc, non-uniform models. 

We thank the ANU supercomputing facility time assignment committee for allocations of computer 
time. This research was also funded by ARC discovery project grants DP0345983 and DPDP0664434. 
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